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Abstract 

Cerebrovascular diseases such as braiu aueurysuis are a priuiary cause of adult disability. The 
flow dyuauiics iu braiu arteries, both duriug periods of rest aud iucreased activity, are kuowu to 
be a uiajor factor iu the risk of aueurysui foruiatiou aud rupture. The precise relatiou is however 
still au opeu field of iuvestigatiou. We preseut au autouiated euseuible siuiulatiou uiethod 
for uiodelliug cerebrovascular blood flow uuder a rauge of flow regiuies. By autouiatically 
coustructiug aud perforuiiug au euseuible of uiultiscale siuiulatious, where we uuidirectioually 
couple a ID solver with a 3D lattice-Boltzuiauu code, we are able to uiodel the blood flow iu 
a patieut artery over a rauge of flow regiuies. We apply the uiethod to a uiodel of a luiddle 
cerebral artery, aud hud that this approach helps us to fiue-tuue our uiodelliug techuiques, aud 
opeus up uew ways to iuvestigate cerebrovascular flow properties. 

Keywords: multiscale modelling, blood flow, ensemble simulation, parallel programming, high-performance 
computing 


1 Introduction 

Stroke is a uiajor cause of death aud luorbidity iu the developed world. Subarachuoid haeuior- 
rhage (SAH) is a type of stroke characterised by bleediug iuto the fluid arouud the braiu, for 
exauiple due to the rupture of au iutracrauial aueurysui. Au aueurysui is a cougeuital weakuess 
iu a blood vessel wall which gradually bulges out to forui a balloou which cau eveutually burst. 
SAHs represeut 5% of cases of stroke, but is relatively luore iuiportaut, as the luortality rate 


1 



Multiscale ensemble blood flow 


Itani et al. 


for these events is about 50%. Overall, approximately 5-10 people per 100,000 are affected by 
SAH due to bleeding in the intracranial arterial wall. [T] The mean age of the victims is 50 years 
and 10-15% fail to reach hospital. Unruptured aneurysms are much more prevalent, estimated 
to affect 1-5% of the population of the UK [2]. Indeed, unruptured / asymptomatic cerebral 
aneurysms are a relatively common finding when scanning the brain for other reasons [T] . Cur¬ 
rent methods of determining which aneurysms have a significant risk of subsequent rupture 
are based on crude measures such as aneurysm size and shape, and there is a clear need for a 
non-invasive tool to stratify risk more effectively in this large patient group. 

Computational fluid dynamics (CFD) techniques may provide means to help quantification of 
the rupture risk, if they can incorporate the key conditions affecting brain aneurysms. Partic¬ 
ularly high or low wall shear stress is believed to increase the risk of aneurysm rupture [3] . Re¬ 
searchers increasingly apply computational fluid dynamics to investigate these problems mum, 
and in particular Shojima et al. concluded that both a very high and a very low wall shear 
stress increases the chance of aneurysm growth and rupture in MCA aneurysms [7]. In align¬ 
ment with these research efforts, we seek to establish computational diagnosis and prediction 
techniques, which may lead to major health benefits and reduce the costs of health care in the 
long term. 

An essential driver for these CFD calculations is the flow solver, and over the last decade 
several sophisticated and scalable solvers have emerged. Within this work we rely on HemeLB 
(described in Section 2.1), which is highly optimized for modelling sparse geometries and has 
unique optimizations which allow it to achieve excellent load balance in the presence of complex 
boundary and in- and outflow conditions [8]. There are several other scalable flow solvers that 
are worth mentioning as well. These include the Nektar finite element package [UllOlin], the 
Palabos package [HI [T3l HI] , the Musubi environment [15] [16] , MuPhy [17] and WaLBerla [18] . 
Although the aforementioned works have provided valuable insight into the haemodynamic 
environment of brain aneurysms, little is known about how the intrinsic variablity of blood flow 
throughout the day affects aneurysm growth and rupture. 

The purpose of this paper is to present a tool which automatically creates an ensemble of 
multiscale blood flow simulations based on a set of clinically measurable patient parameters, 
and runs these simulations using supercomputing resources. The tool allows us to automate 
the study of the blood flow in a vascular geometry under varying patient-specific conditions. 
In addition, an automated data processing component extracts velocity and wall shear stress 
(WSS) values, and generates plots and animations which allow us to visualize these properties 
in the vascular geometry. This paper builds on previous works where we simulated flow in 
arterial networks using a single flow configuration [ISEolET]. 

To showcase our approach, we construct and execute a range of multiscale simulations of a 
middle cerebral artery (MCA). We present the results of these simulations, and compare our 
approach to related efforts as an initial validation of our 1D-3D multiscale scheme. This work is 
organised as follows. In Section 2, we present the tools we developed to perform our multiscale 
ensemble simulations and how we integrate them in an automated workflow. We describe the 
setup of our simulation in Section 3, our results in Section 4 and provide a brief discussion in 
Section 5. 


2 Automated multiscale ensemble simulations 

Our automated workflow combines three existing components. These include the HemeLB and 
pyNS simulation environments, and the FabHemeLB automation environment. In this section 
we describe these three components, and how they interoperate in our automated multiscale 
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ensemble simulation (MES) environment. 

2.1 HemeLB 

HemeLB is a 3 dimensional lattice-Boltzmann simulation environment developed to simulate 
fluid flow in complex systems. It is a MPI parallelised C++ code with world-class scalability 
for sparse geometries. It can efficiently model flows in sparse cerebral arteries using up to 
32,768 cores [22l [23] and utilises a weighted domain decomposition approach to minimize the 
overhead introduced by compute-intensive boundary and in-/outflow conditions [8|. HemeLB 
allows users to obtain key flow properties such as velocity, pressure and wall shear at predefined 
intervals of time, using a property-extraction framework. 

HemeLB has previously been applied to simulate blood flow in healthy brain vasculature as 
well as in the presence of brain aneurysms pnipi]. Segmented angiographic data from patients 
is read in by the HemeLB Setup Tool, which allows the user to visually indicate the geometric 
domain to be simulated. The geometry is then discretized into a regular unstructured grid, 
which is used as the simulation domain for HemeLB. HemeLB supports predefined velocity 
profiles at the inlets of the simulation domain, which we generate using pyNS in this work. 

2.2 pyNS: Python Network Solver 

pyNS is a discontinuous Galerkin solver developed in Python, which simulates haemodynamic 
behaviour in vascular networks [25]. pyNS uses aortic blood flow input based on a set of patient- 
specific parameters, and combines one-dimensional wave propagation elements to model arterial 
vasculature with zero-dimensional resistance elements to model veins. The solver requires two 
XML files as input data, one with a definition of the vasculature and one containing the simu¬ 
lation parameters. Simulation parameters include mean blood pressure, cardiac output, blood 
dynamic viscosity and heart rate. pyNS has been used in several studies, e.g. to try to inform 
treatment decisions on haemodialysis patients [26] and as a large-scale model for distributed 
multiscale simulations of cerebral arteries Hg. 

2.3 FabHemeLB 

FabHemeLB is a Python tool which helps automate the construction and management of en¬ 
semble simulation workflows. FabHemeLB is an extended version of FabSim [27] configured 
to handle HemeLB operations. Both FabSim and FabHemeLB help to automate application 
deployment, execution and data analysis on remote resources. FabHemeLB can be used to com¬ 
pile and build HemeLB on any remote resource, to reuse machine-specific configurations, and 
to organize and curate simulation data. It can also submit HemeLB jobs to a remote resource 
specifying the number of cores and the wall clock time limit for completing a simulation. The 
tool is also able to monitor the queue status on remote resources, fetch results of completed 
jobs, and can conveniently combine functionalities into single one-line commands. In general, 
the FabHemeLB commands have the following structure: 

fab <target machine> <command>:<parameter>=<value>,... 

For example: 

fab archer ensemble:config=/path/to/config,cores=1536,wall_time=05:00:00 

In table we present a number of commands typically executed with FabHemeLB in the scope 
of this work. The commands are customised to run on local machines, continuous integration 
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Table 1: List of commands commonly used in FabHemeLB. 


command name 

brief description 

cold 

Copy HemeLB source to remote resource, compile and build every¬ 
thing. 

run_pyNS 

Execute instances of pyNS to generate a range of flow output files. 

generate_LB 

Convert pyNS output to HemeLB input. 

submit _LB 

Given a set of velocity profiles, submits the corresponding HemeLB 
jobs to the remote (supercomputer) resource. 

fetch_results 

Eetch all the simulation results from the remote resource and save 
them locally. 

analyze 

Performs data-analysis that allows for easy visualization of the re¬ 
sults. 

ensemble 

Do all of the above, except cold. 



Figure 1: Workflow diagram showing the processes involved in the ensemble simulation method. 
The simulations were distributed, with PyNS simulations using a local workstation in London, 
and the HemeLB simulations using the ARCHER supercomputer. 


servers, regional, national or international supercomputing resources. The workflow is presented 
in the diagram of Figure Specifically, the analyze command processes the compressed output 
files to generate human readable files and visualisations. It also generates an image file showing 
the whole geometry, wall shear stress within the geometry, and velocity measurements on pre¬ 
selected planes inside the geometry over time as an animation. 
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3 Setup 

To apply our automated ensemble simulation tool, we employ patient-specific parameters mea¬ 
sured by Sugawara et al. during a study to assess cardiac output during exercise [28]. They 
measured the blood pressure, cardiac output and heart rate of 16 young patients (9 male and 
7 female) at different exercise intensities, being at rest or at 70%, 90%, 110% and 130% of the 
ventilatory threshold (VT). The VT is the point during exercise training at which pulmonary 
ventilation becomes disproportionately high with respect to oxygen consumption. This is be¬ 
lieved to reflect onset of anaerobiosis and lactate accumulation. We add two more sets of values 
at 30% and 50% VT, linearly extrapolating the other parameters. The resulting seven sets of 
parameters are used in our automated workflow, resulting in seven HemeLB simulations being 
run. We present the seven sets of parameters in Table 


Table 2: Configuration data used as input for pyNS to run the ensemble of simulations. The 
values are based on the average of the measurements of 16 people at different exercise intensities 
measured by the percentage of the ventilatory threshold(VT) [28]. In the last column, we pro¬ 
vide the mean flow velocity in the right MCA, as calculated using PyNS, for each configuration. 


Configuration 

Exercise 

Blood Pressure 

Cardiac 

Heart 

mean flow 


intensity 

Mean 

output 

rate 

velocity 



[mmHg] 

[L / min] 

[bpm] 

[ms~^] 

1 

Rest 

80 

4.8 

68 

0.460 

2 

30% VT 

87 

6.2 

79 

0.451 

3 

50% VT 

94 

7.6 

90 

0.428 

4 

70% VT 

100 

9 

101 

0.393 

5 

90% VT 

112 

10.7 

113 

0.371 

6 

110% VT 

116 

11.9 

120 

0.351 

7 

130% VT 

122 

13.2 

134 

0.339 


For our pyNS simulations, which we ran for 10 cardiac cycles, we set the blood density to 1050 
Kgjnn?^ Poisson’s ratio of transverse to axial strain to 0.5 and the time step to 5 ms. For 
our HemeLB runs we use a model derived from a patient-specific angiographic 3D geometry of 
a middle cerebral artery, supplied by the Lysholm Department of Neuroradiology, University 
College Hospital, London and segmented using the GIMIAS tool [29]. We use a voxel size of 
18.9/im, which results in a geometry containing 13,179,961 lattice sites. In Figure we show 
the setup tool interface with the MCA geometry. Here the inlet is given by the green plane and 
the outlets by the red planes. The location of interest for our WSS analysis is highlighted. For 
simulations of this particular voxelized simulation domain, we specify a time-step of 0.5014/i5 
and run each simulation for 7.9 million steps, or 4 seconds of simulated time. 

During our HemeLB runs, we store the WSS throughout the geometry for every 50,000 time 
steps. In addition, we define an output plane close to the outlet, at 49mm from the ear, 
2mm away from the outflow boundary (or outlet)^ to record velocity and pressure data every 
50,000 steps. In all our runs we use interpolated Bouzidi wall boundary conditions m and 
zero-pressure outlet conditions, (for details see [22]). Using the ensemble command, we run an 
ensemble of 7 simulations on the ARCHER supercomputer. In total, we used 10,752 (7 x 1536) 
cores for a duration of approximately 3.5 hours. 
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Figure 2: HemeLB Setup Tool loaded with a middle cerebral artery geometry. The inlet is 
indicated by a green surface, and the three outlets with a red surface. The lattice site which 
we used for in-depth wall-shear stress analysis is highlighted with a light blue arrow. 


4 Results 

We present the time series of the maximum velocity in the measurement plane in Figures 
and for the parameters listed in Table The curves demonstrate that the frequency of 
the ID cardiac cycles generated by pyNS are accurately reproduced in the 3D high-resolution 
HemeLB simulations. The peak velocity in the measurement plane is higher than the input 
velocity because the measurement plane has a lower area but the total flux remains constant 
throughout the MCA due to the incompressible flow. 

Within the ensemble, the frequency of the cardiac cycles increases with exercise intensity as the 
heart rate increases. While this seems trivial, it indicates that the time-resolution chosen to 
couple pyNS and HemeLB is sufficient to reproduce this effect. The good match between the 
input and the measurement gives confidence that the LB simulations have converged and the 
cardiac cycles are stable. 



Figure 3: Velocity time-series for the configurations 1 (bottom lines) to 3 (top lines). Dashed 
lines indicate the velocity time-series generated for the inflow boundary by pyNS, and solid 
lines correspond to the maximum velocity measured at the 49mm plane in HemeLB. 
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Figure 4: Left: Velocity time-series as in Fig. but for configurations 4 (bottom line) and 5 
(top line). Right: Velocity time-series as in FigTj^but for configurations 6 (bottom line) and 7 
(top line). 


4.1 Wall shear stress 

Figure shows snapshots of the WSS for four configurations of the ensemble at the peak sys¬ 
tole, when the flow velocity is highest. The colour scale allows us to identify regions of high 
WSS which are located predominantly at constrictions of the MCA and around the inlet in 
an assymetric pattern. Based on the WSS values near the inlet, we conclude that the stan¬ 
dard circular-shaped Poisseuille velocity profile used in HemeLB is ill suited for patient-specific 
geometries, as they usually feature non-circular inlets. As a result, we are now developing a 
method for a modified velocity profile which takes non-circular inlet shapes into account. 

At the point of interest indicated in Fig. we observe much higher wall shear stress at higher 
exercise intensity. We present a cross-instance analysis of the WSS at this point in Fig.[^ Here 
we provide for instance the mean WSS, which decreases linearly with the mean velocity at the 
inlet (see Tablej^ both values are reduced to ^0.75 times the magnitude at rest, when measured 
at 130% VT exercise intensity). This matches the theoretically expected linear scaling of the 
WSS with the velocity parallel to the wall. This parallel velocity is expected to be proportional 
to the velocity at the inlet for simple geometries and flow regimes. Our WSS results are in line 
with related literature, which report maximum values in MCAs in the range of 14 to 40 Nm~‘^ 
for unruptured aneurysm geometries mm- 

We also present the maximum WSS in time, which is ^ lSNm~^ for the run at rest, with 
a heart rate of 68 bpm and a maximum velocity of 0.84ms“^. At full exercise intensity, the 
maximum WSS is much higher, at ^ 31 Nm~^ for a heart rate of 134 bpm and a maximum 
velocity of 1.19m5“^. Here, while the mean WSS and velocity decrease with exercise intensity, 
the maximum WSS and velocity increase. Between these two cases, the difference in maximum 
WSS (a factor of 1.77) cannot be justified solely by the difference in maximum velocity (a factor 
of 1.45). Indeed, the (much larger) difference in heart rate (a factor of 1.97) may be an important 
contributor to the magnitude of the maximum WSS. Further investigations are required to 
explore the exact nature of these relations. The variability of the WSS, here measured as the 
averaged absolute difference between consecutive extractions at a 0.025 s interval, increases as 
expected with higher exercise intensity. 
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Figure 5: A snapshot of the wall shear stress for the configurations 1, 3, 5 and 7 at the peak 
systole. Red regions indicate high wall shear stress while blue indicates low wall shear stress. 
Exercise intensity is increasing clockwise for the configurations shown. 



Figure 6: Cross-instance analysis of the WSS for the location of interest indicated in Fig.[^ We 
extracted WSS values in the simulations for a range of 3 cardiac cycles (when at rest) up to 6 
cardiac cycles (at 130% VT). We present the maximum measured WSS (blue line), the average 
WSS (green), and the time average of the WSS slope, calculated over intervals of 0.025 s (red). 
The WSS values at rest can be found on the left side (0% VT). 


5 Discussion and conclusions 

We present an automated ensemble simulation framework and its application to model blood 
flow in the middle cerebral artery under a range of patient-specific cardiac parameters, using a 
multiscale ensemble approach. We show good agreement of velocity profiles at the inlet with 
those close to the outlet, and that our non-lattice aligned inflow conditions require further 
enhancement. FabHemeLB allows us to run the whole workflow for the relatively complicated 
setup in one tool, including the execution and analysis of the ensemble simulations. It reduces 
the human effort required for doing these tasks, and by automatically scheduling the ensemble 
instances in parallel it also allows for efficient use of large core counts and a reduced time to 
solution. The systematic execution and analysis patterns offered by FabHemeLB allow us to 
easily identify shortcomings in our existing approach. Not only does this feature in FabHemeLB 
boost our ongoing research, it also provides the level of data curation required to do future. 
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more extensive, validation studies. 

In our case study, we investigate the wall shear stress (WSS) properties in a middle cerebral 
artery at a location of interest close to the outlet. We find that the mean WSS correlates 
as expected linearly with the average flow velocity at the inlet. However, in addition we find 
evidence that the maximum WSS is dependent on the heart rate as well as the average flow 
velocity. This implies that these relations are non-trivial, and that a comprehensive analysis 
of flow dynamics in cerebral arteries should not only include the presence of pulsatile flow, but 
also the presence of these flows over a range of heart rates. 
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